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ALGORITHMS  FOR  COMPUTATION  OF  WATER  LEVEL  ELEVATION 
AT  A  FIXED  LOCATION  FROM  THE  WATER  LEVEL 
ELEVATIONS  AT  A  MOVING  PLATFORM 


by  Leon  Borgman 


INTRODUCTION 

The  basic  input  data  for  this  problem  consists  of  water  level  elevations  measured  at  a 
succession  of  times  from  a  moving  platform.  Let  tf,  t2,  t3,  •  •  •,  tj^  be  the  times  at  which 
measurements  are  made.  Positioning  devices  are  used  to  establish  the  platform  positions 
(x(tn),  y(tn);  n  =  1,  2,  •  •  •,  M}.  The  corresponding  sequence  of  water  level  elevations  at 
the  platform  are  {7^,  n  *  1,  2,  •  •  •,  M}.  For  the  reference  position  (xo,  yo),  let  7o(tn)  be 
the  water  level  elevation  at  the  reference  location  at  time  tn.  The  basic  problem  is  to 
predict  {j?o(tn),  n  «  Mi,  Mj+1,  •  •  •,  M2}  for  1  <  Mi  <  M2  <  M  from  {ijn,  n  =  1,  2.  •  •  •, 
M}. 

This  problem  can  be  investigated  with  geostatistical  techniques  if  the  directional 
wave  spectrum  S(f,  ff)  is  available  from  nearby  measurements.  Here,  and  in  the  above, 

f  =  frequency  in  cycles  per  second,  Hertz. 

x  as  a  horizontal  direction. 

y  —  a  second  horizontal  direction  perpendicular  to  the  x— axis. 

9  =  an  angle  measured  from  the^  positive  x— direction  toward  the  positive 
y— direction  so  that  the  positive  y— axis  is  at  9  =  90°.  (Note:  this 

definition  of  9  is  (x,  y)  —  axes  dependent,  may  differ  markedly  from 
compass  directions,  and  may  end  up  clockwise  or  counter-clockwise.) 


2 


S(f,  8)  =  a  spectral  density  function  defined  so  that 

rto  (*2t 

2  S(f,0)dfldf  =  <r2 

J  A  J  A 


0^0 


where 


<r2  =  variance  of  sea  surface  elevation. 

Algorithms  for  estimating  rj o(t)  by  the  two  geostatistical  techniques  of  simple  linear 
kriging  and  conditional  simulation  will  be  developed  in  the  following  report.  Both  of  these 
techniques  depend  on  the  covariance  matrix  of  the  vector 

71  1 

72 


% 

70>M, 

7o»M2 

Hence,  some  time  will  be  spent  initially  defining  the  covariance  matrix  and  developing 
procedures  by  which  it  can  be  computed  from  S(f,0).  This  will  be  followed  by  two  sections: 
the  first  on  kriging  methods,  and  the  second  on  conditional  simualtion.  Finally  a  summary 
and  conclusions  section  will  be  given. 


» 

a 

.  Ho  . 
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The  algorithms  for  the  kriging  and  conditional  simulation  are  quite  straightforward 
once  the  covariance  matrix  is  developed.  The  main  task  to  be  completed  is  the 
development  of  a  rapid  procedure  for  computing  the  covariances.  The  major  part  of  the 
report  will  be  concerned  with  this  problem. 


Let  (xt,  Ji)  and  (X2,  72)  be  an7  two  locations  and  (tj,  to)  be  an7  two  times.  The 
water  level  elevation  above  mean  water  level  at  (x,  7,  t)  will  be  denoted  b7  t](x,  7,  t). 

The  covariance  between  two  random  variables  U  and  V  is  defined  as  the  statistical 
average  of  (U  —  fMi)(V  —  /*v)  where  /i<,  is  the  statistical  average  of  U  and  y-,  is  the  statistical 
average  of  V.  It  is  convenient  to  denote  the  statistical  averaging  process  with  E  [  ].  Then 
the  covariance  definition  ma7  be  expressed  in  S7mbols  as 

Cov(U,  V)«E[(U-/*XV-*)] 

where 

B(  U  ]  —  Ma 
E IV]-* 

Now  replacing  U  with  ij(x  1,  j\,  t()  and  V  with  r](x 2,  72,  t2),  and  noting  that,  b7  definition 
of  17,  the  average  of  jj(xi,  7i,  t|)  and  tj(x 2,  72,  t2)  is  zero,  the  covariance  between  the  two 
can  be  written  as 

Cov[jj(xi,  71.  t|),  Ij(x2,  72,  t2)]  =  Efi^xt,  71,  t|)  fix 2,  72,  t2)l 

That  is,  the  covariance  is  the  expected  product. 

It  will  be  asusmed  that  the  sea  surface  is  a  stationar7,  Gaussian  process.  Then  the 
covariance  defined  above  will  onl7  be  a  function  of  the  differences  in  position  and  time.  Let 


T  =  t2  —  t| 

X  =  *2  —  X| 

Y  =  y2  -yi 

be  these  differences,  and  define  the  covariance  function  as  C(X,  Y,  r)  where 


C(X,  Y,  r)  =  Cov[ij(xi,  yi,  t,)(  jXx2,  72,  t2)] 


It  can  be  shown  that  (Borgman,  1969,  p.  723)  the  spatial,  temporal 
function  can  be  expressed  in  terms  of  the  directional  spectrum  as 


C(X,  Y, 


S(f,0)  cos[wk(Xcos0  +  Ysm0)  —  2rfr]d0df 


where 


w 


+1,  if  0  is  the  direction  toward  which  waves  travel 
— 1,  if  9  is  the  direction  from  wh  i  ch  waves  travel 


k  =  wave  number 

2 

»  function  of  f  defined  by  (2xf)  as  gk  tank  (kd) 
d  as  water  depth 
g  =a  acceleration  due  to  gravity. 


The  functions  involved  can  be  defined  for  negative  frequency  by 
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covariance 


S( — C,  ff)  =  S(f,  &) 
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k(-f)  =  -4(f) 


Then,  through  the  complex— valued  definition  of  the  cosine  as 


cos  <j>  =  [exp(i  <f>)  —  exp(— i  d)]/2 


the  formula  for  C(X,  Y,  r)  can  be  rewritten  as 


C(X,  Y,  r) 


S(f,  0)  exp[— iwk(Xcos0  +  Ysin0)  +  i2xfr]d#df 


Suppose  (X,  Y)  is  re— expressed  in  polar  coordinates  ( p ,  a),  with 


X  =  p  cos  a 


Y  =  p  sin  a 


then 


X  cos  9  +  Y  sin  9  =  p  cos  a  cos  9  +  p  sin  a  sin  9  =  p  cos  ( 9— a) 


With  these  definitions 


C(X,  Y 


,  r)=r  rs 

J  J  A 


S(f,0)  exp[— iwk  /?cos(d-or)]  exp[i2*fT]d&if 


Without  loss  of  generality,  S(f,0)  can  be  expressed  in  the  product  form 


S(f,  0)  =  S(f)Df(0) 


where  Df(  9)  is  the  spreading  function  defined  so  that 


r2x 

Df(^)d0=  1.0 
J  o 
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The  spreading  function  gives  the  distribution  of  wave  energy  or  variance  with  direction  at 
frequency  f.  The  covariance  function  can  be  written  in  terms  of  the  product  form  as 


C(X,  Y,  r)  =  J°°  S(f)  |  ^  Df(  &)  cos  ( ^dfl} 


*ei2jfrdf=C*(/>,a) 


(where  *  denotes  multiplication). 

In  summary,  the  covariance  can  be  expressed  in  rectangular  and  polar  form  as 


C(X,  Y,  r)  =  2  f*1  S(f)  |  f  Df(0)  cos[wk(Xcos0  +  Ysintf)  —  2Tfr]d^|df 
o  1  -1  o  * 

C *(/>,  a,  r)  =  J"  S(f)  {  J2'  D,(6)  e-i*k,w»(«-a)d,j  ei2rfrd( 

Either  form  can  be  used  as  a  basis  for  computaitons.  In  practice,  it  is  usually  best  to  pre 
compute  a  table  of  the  Covariance  function  values  for  a  given  S(f,  9)  and  a  grid  of  values, 

—ho  £  X  <  ho 
— ho  <  Y  <  ho 

0  <  T  <  To 

or 

0  <  p  <  po 
0  <  a  <  2t 
0  <  r  <  to 


It  is  only  necessary  to  compute  these  for  positive  time  lag  because,  from  the  definitions 


C(— X,  -Y,  -r)  =  C(X,  Y,  r) 

C*(p,  a+T,  -r)  =  C*(p,  a,  r) 

For  the  Navy  application,  where  the  S(f,0)  function  is  estimated  from  a  buoy,  further 
simplications  are  often  appropriate.  It  b  often  satisfactory  to  take  Df(  9)  as  only  depending 
on  9.  When  this  happens,  the  spreading  function  will  be  denoted  by  D(  0).  Three  common 
formulas  for  the  spreading  function  that  are  used  as  approximations  are 
The  generalized  cosine— squared  model, 

D(0)  =  ccos2*[(0-*i)/2] 


The  von  Mises  model, 


D(0)  =  (exp(acos(^-/i)]}/{2Tlo(a)} 


where  Io(a)  =  modified  Bessel  function  of  order  zero. 


and  the  wrapped— normal  model 


D(0)  =  J  exp 


]  =  -T30 


/2 


Kjn  a) 


All  models  have  about  the  same  shape  and  are  unimodai  and  symmetric  about  /± 
Approximate  equivalent  values  between  s,  a,  and  a2  are  given  in  the  appendix  to  the 
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report. 

Since  all  three  models  have  very  nearly  the  same  shape  for  most  typical  wave 
spreading  conditions,  it  is  just  a  matter  of  mathematical  convenience  which  is  used.  If  one 
is  a  reasonable  approximation,  then  any  other  of  the  three  will  also  be  a  reasonable 
approximation. 

In  the  polar  covariance  function  with  the  spreading  function  independent  of 
frequency,  it  is  particularly  convenient  to  use  the  von  Mises  model. 


C  *(/>,<*,  r)  = 


_r 


J. 


s(o{  r 


■2x  gacos(0-^)  -  iwkpcos(#-or) 


dd 


2xIo(a) 


eErfrdf 


There  is  a  nice  series  approximation  of 


eacos(  $-p) 

in  terms  of  the  modified  Bessel  functions  given  by  Oliver  (1964,  p.  376,  eq.  9.6.34)  as 

eacos (9-ii)  =  I(j(a)  +  2  |  ln(a)cos[n(  £-/*)] 

n=l 

After  a  little  algebra,  the  formula  for  C *  can  be  expressed  in  series  form  as 
C*(/>,  «,  r)  »  k  r  S(()  (a'e-<wW*<»(«-a)djei2rfr<lf 

■'ll)  ■*  o 

+  S(f)  J2T  cos(n(^)]e-iwk^os(  ^W^df 

The  integration  over  (0,  2x)  is  really  just  an  integration  over  the  full  circle  of  360°. 
The  full  circle  integration  could  just  as  well  range  over  (a  —  r,  a  +  r).  If  the  limits  of 
integration  are  changed  to  this  new  choice  and  the  variable  of  integration  is  changed  to 

i>  =  9  —  a 

the  formula  simplifies  to 

C‘(p,  a,  r)  =  f  S(f)  f%^wk^oa  W2rfrdf 


+  I  j”  s(0  f  CO ,i2rfrd{ 

The  cosine  in  the  second  expression  can  be  expanded  to 

cos[n(^— ^+Qf)]  =  cos(n  y)cos[n(  o— ^)]  —  sin(ni^)sin[n(  a—fi)} 

With  this, 

C*(*  a,  r)  =  L  f  S(f)  f 

+  S  C03(n(  <*-^1  s(f)  f  cos(n^) 

n=l  v  '  J-oo  t 

*e-4wkpcos^d^ei2xfrdf__  J  sin[n( «-*)J  j*  S(f)  sm(ntf) 

*e-iwk/x:os^d^ei2T{Tdf 

The  last  expression  is  zero  since  the  integrand 

3in(n*.) 

is  an  odd  function  of  rp  integrated  over  (—r,  r).  The  othertwo  integrands 

g— iwk/xos^ 
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and 

cos (n*) 

are  even  functions  of  i>,  so  the  range  of  integration  can  be  changed  to  (0,  t)  with  a 
multiplication  by  2.  Thus 

C*U  a,  r)  =  I  f  S(f)  [T  e-iwk^°^d^  eiMrdf 

J-ffl  J  Q 

S(t)JTcos(ad') 

,e-!wkpcos^ei2*{rd{ 

The  reason  for  all  these  manipulations  will  now  become  apparent,  the  integrals  over 
i)  can  be  expressed  as  Bessel  functions.  Oliver  (op.cit.,  p.  360,  eq.  9.1.21)  gives 

f  g-iwk/xos^^n^^  _  inTJn(— wk/?) 

Jo 

r' ^wkpcos*^  = 

J  0 

By  Oliver  (op.cit.,  P.  360,  eq.  9.1.20)  if  the  argument,  z,  is  real— valued 


JwM  =  (-DQ  Jn(») 
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Combining  all  these  results 


C*(p,  a,  r)  =  f  S(f)  Hip)  e'2rirdf  +  2  S  (-»>)“  Mil  cos[n(o-,i)l 

n=l  l0^a' 

•f  S(f)  J„(Mei2rfTdf 


The  last  equation  is  in  a  form  appropriate  for  approximation  with  the  fast  Fourier 
Transform.  Let  N  and  Af  be  chosen  so  that  S(f)  Jn(k^)  is  essentially  zero  for  f  >  NAf/2 
and  n  =  0,  1,  2,  3,  4,  •••.  This  is  equivalent  to  choosing  a  frequency  beyond  which  S(f) 
will  be  treated  as  though  it  is  exactly  zero.  Then  define 

Ars(NAf)’1 

A(a)  =  S(mAf)  Jn(kn  p) 

la 

where  km  is  the  wave  number  corresponding  to  f  =  mAf,  and  set 


4°J>=a(:> 


Then,  introducing  a  new  function,  R 0(p,  r ) 


Rn(p,  jAr)  =  r  S(f)  JKk/?)ei2rf^Ar)df  s  Af  A(a)el2^m/N 
J-"'  m=0 


the  function  Rn(/>,  r)  can  be  computed  quite  rapidly  for  a  selected  list  of  p— values  to 
develop  a  matrix  whose  rows  are  the  p— values  and  whose  columns  are  the  r  =  jAr 
time-lag  values. 

An  algorithm  for  the  rapid  computation  of  Jn(z)  for  real— valued  z  is  given  by  Oliver 
(op.cit.,  bottom  of  page  385).  An  exactly  parallel  procedure  can  be  applied  to  compute 
In(a)  with  eq.  9.6.36  (Oliver,  op.cit.) 
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In  terms  of  the  function  Rn(/?,r)  defined  above 

C*u  O,  T)  =  RoU  r)  +  2  !  (-wi)°  A  r)  cos[n(  a-).)] 

n=l  i<na' 

The  nature  of  values  of  Rn(yt>,  t)  can  be  combined  with  any  a  to  compute  rapidly  the  value 
of  a ,  r). 


Let  Cu  be  the  M  by  M  matrix  whose  (nj,  112)  element  is  the  covariance  between  rjn^ 
and  rjn^  for  the  water  level  elevations  in  the  series 


{%,  n  =  1,  2,  •  •  •,  M} 


defined  in  the  introduction.  Similarly  let  C22  be  the  (M2  —  Mi+1)  by  (M2  —  Mj+1)  matrix 

whose  (ji,  \t)  element  is  the  covariance  between  and  77o( t j -|- j 2 1 ) ’  ‘n 

series  {qo(tn),  n  =  Mi,  M|+l,  •  •  •,  M2}  defined  in  the  introduction.  Finally  let  C12  be  the 
M  by  (M2  —  Mi+1)  matrix  whose  (n,  j)  element  is  the  covariance  between  7^  and 
’^Mi+j— these  definitions,  the  covaraince  matrix  of  (tj,  2Z0)  is 


Cov 


‘  n  ' 

C  ,  1  C,2  * 

— ■ 

.  CT2  C22 . 

All  of  the  covariances  in  Cu,  C12,  and  C22  can  be  computed  by  first  determining  X, 
Y,  and  r  for  the  pair  of  locations  and  times  involved,  then  getting  p  and  a  from 

/>  =  /  X7  +"  Y7 


a  =  arc  tan  (Y/X), 


and  finally  interpolating  for  R n(p,  r)  and  computing  C*(p ,  a,  r). 
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CONDITIONAL  SIMULATION  OF  ^ 

The  matrix  formula  for  conditional  simulation  of  given  the  values  of  rj,  was 
derived  previously  by  Borgman  (1984,  p.  533,  eq.  85.)-  The  following  steps  are  involved. 


(1)  Compute  the  eigenvectors  and  eigenvalues  of 


C  = 


C 1 1  C  12 
C  |2  C22 


Let  Aj  be  the  eigenvalues  and  Vj  be  the  corresponding  eigenvectors.  Define  V  as  a  matrix 
whose  columns  are  the  eigenvectors  and  L  as  the  diagonal  matrix  whose  main  diagonal 
elements  are  Aj  and  whose  off-diagonal  elements  are  zero.  Then 


C  =  VLVT 


in  the  well  known  eigenvector,  eigenvalue  decomposition  (Jennings,  1977,  p.  32,  eq.  1.130). 
Consequently, 


C  =  (VL'/2)(VL,/2)T 


(2)  Let  2l  be  a  vector  of  independent,  standard  normal  random  numbers.  Then 


it 


Ho 


=  VL1/2 Z* 


is  an  unconditional  simulation  of  3  and  rj^. 
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(3)  A  conditional  simulation  of  r^,  given  the  values  of  ij  is 


{.2o>  n)  = 


c!,  c;i(a-3*)  +  n, 
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KRIGING  FORMULAS  FOR 

The  kriging  procedure  estimates  each  7o(tj)  by  a  linear  combination  of  the  interval  of 
values  of  ffo  surrounding  the  time  tj.  The  idea  here  is  to  use  all  values  of  that  are 
correlated  with  7o(tj).  A  reasonable  choice  would  be  to  include  all  7q  measured  from  the 
platform  within  two  wave  periods  of  the  time  tj.  Here  a  good  value  for  the  wave  period 
would  be  the  inverse  of  the  frequency  at  the  peak  of  the  spectra.  The  value  rjo(tj)  is 
estimated  by 

%(*\)  =  E  anrjn 

where  the  summation  extends  over  the  times  surrounding  tj.  The  coefficients,  an,  are 
computed  to  minimize 

Q  =  E[{t?0(tj)-  Tk(tj)}2] 

subject  to  the  constraint  that 

E(»?0(tj)-»?0(tj)]  =  0 

A  full  discussion  of  this  procedure  i3  given  by  Bergman  (1985,  p.  14),  a  copy  of  which  is 
forwarded  with  this  report. 

The  computations  may  be  summarized  as  follows.  Let  ni  <  n  <  n?  be  the  interval  of 
values  of  7 n  used.  » 
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C 1 1  =  (no— ni+1)  by  (112— fli+l)  covariance  matrix  of  the  ^  for  ni  <  n  <  no 


Cm  =  vector  whose  f— th  element  is  the  covariance  between  ^  and  70(tj) 
^  =  vector  of  n2— m+1  components,  all  of  which  equal  1.0 


a  =  vector  of  coefficients  to  be  multiplied  by  the  tj n  in  the  interval 


The  kriging  equations  which  solve  the  constrained  minimization  are  computed  from 


Cm  1 

a " 

C 1  n 

aT  oj 

.  A. 

1.0 

where,  A  is  the  Lagrangean  multiplier  imposing  the  constraint.  Thus 


’  a 

■ct, 

1 ' 

-1 

A 

,iT 

0. 

.1.0. 

The  mean— square— error  of  the  estimate  of  70(tj)  is 

mean  square  error  =  E[{^0(tj)  —  70( t j  )}2]  =  a2  —  A  —  a^Cm 

where 


<72  =  2  T  S(f)df 
J 


from  the  sea  surface  spectral  density. 
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SUMMARY  AND  CONCLUSIONS 

1.  Kriging  and  conditional  simulation  algorithms  have  been  developed  for  estimating 
the  water  level  elevation  at  a  fixed  reference  location  from  the  measured  water  level 
locations  on  the  moving  platform. 

2.  Both  kriging  and  conditional  simulation  have  straightforward  mathematical 
formulas,  once  the  appropriate  covariance  matrices  have  been  computed. 

3.  An  algorithm  based  on  a  Bessel  function  series  and  the  use  of  the  fast  Fourier 
Transform  to  compute  a  subsidiary  function,  Rn(p,  r),  is  derived.  This  is  the  main  work  in 
the  report.  Both  kriging,  and  conditional  simulation  procedures  are  completely  derived  in 
cited  references. 

4.  Recommendations  are  given  for  procedures  to  compute  the  Bessel  functions  in  the 
formula  and  the  other  aspects  of  the  algorithm. 

5.  The  kriging  procedure  will  probably  be  the  best  choice  of  the  geostatistical 
techniques  for  use  in  estimating  the  water  level  elevations  at  the  reference  location.  A 
convenient  error  measure  arises  naturally  from  the  computations  and  the  procedure  is 
relatively  rapid  to  compute. 
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APPENDIX 

The  equivalences  between  the  spreading  function  models  are  most  easily  derived  from 
the  half— peak  width  of  the  functions.  The  half— peak  with  of  D1&)  is  two  times  the  (B—n) 
value  at  which 

D(0*)  =  O.5DO i) 

Let  Ajjp  be  this  value 

^HP  =  ^ 

For  the  generalized  cosine— squared  model 

D(/i)  =  c 

so  the  equation  becomes 

c  cos2»[(0*— fi)(2]  =  0.5c 
cos[(0*-/i)/2]  =(0.5)‘/(2s) 

A^p  =  2(  9+—fi)  =  4  arc  cos[0.5  ^2s^] 

For  the  von  Mises  model 


D(/i)  =  I/{2xI0(a)} 


25 


Hence  the  equation  is 

eacos(  _  q>5 

acos(  B+— fi)  =  Log(0.5) 

A™  =  2(  ^-/i)  =  2  arc  cos{  ^ M  J 

Finally,  for  the  wrapped— normal  model,  in  the  case  most  common  in  ocean  wave 
work  where 


m 


e-(MV(2»*) 

V5x  <r 


D(/0- 


1 

/Jr  a 


and  the  equation  to  be  solved  is 


e-<  *-#•)’/( 2"2)  =  0.5 


-«^)2 

2a2 


Log(0.5) 


(9+-/i)2  =  — 2<r2Log(0.5) 


Agg  «  2(0*-/i)  =  2/3OTog(03T 
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Thus,  a  reasonable  equivalence  between  parameters  is  given  by 


2  arc  cos 
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